import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import griddata


def test_calculate_mutual_inductance(x1, y1, x2, y2, gapx):
    """
    used to get machine parameters from machine.xlsx file
    path ='C:\\Users\\xiang\\Documents\\XINAO\\AI\\SEPY\\data\\EXL50U\\machine.xlsx'

    input:
        option name         data_type    size    unit    meaning
        1      x1           float        1               x-coordinate of a point in the coil L1
        2      y1           float        1               y-coordinate of a point in the coil L1
        3      x2           float        1               x-coordinate of a point in the coil L2
        4      y2           float        1               y-coordinate of a point in the coil L2
        5      gapx


    output:
        option  name                data_type   size    unit    meaning
        1       mutual_inductance   folat       1       H       mutual inductance between two current elements

    reference:
        [1]皇甫国庆.两圆线圈间互感及耦合系数讨论[J].渭南师范学院学报,2015:26-31.

    Developed by xiangguchn 2023/11/28
    If you have some question, don't hesitate to ask xiangguchn@163.com for help.
    """
    ran_x = np.linspace(0, 60, 60)
    ran_x_1 = np.linspace(0, 20, 20)
    X, Y = np.meshgrid(ran_x, ran_x_1)
    mi = 0 * X
    for i in range(len(ran_x)):
        for j in range(len(ran_x_1)):
            x = ran_x[i]
            x_1 = ran_x_1[j]
            x1 = x2 * x
            y1 = x2 * x_1 + y1
            R = np.sqrt((1 + x) ** 2 + x_1**2)
            b = np.sqrt(4 * x / ((1 + x) ** 2 + x_1**2))
            cmu = np.pi * 4.0e-7
            myk = special.ellipk(b**2)
            mye = special.ellipe(b**2)
            mutual_inductance = cmu * x2 * R * ((1 - b**2 / 2) * myk - mye) * 25

            mi[j][i] = mutual_inductance

    fig = plt.figure("3D Surface", facecolor="lightgray")
    plt.title("3D Surface", fontsize=18)
    ax3d = fig.add_subplot(111, projection="3d")
    ax3d.set_xlabel("x")
    ax3d.set_ylabel("x_1")
    ax3d.set_zlabel("M")
    plt.tick_params(labelsize=10)
    ax3d.plot_surface(X, Y, mi, cmap="jet")
    plt.show()

# mi = test_calculate_mutual_inductance(x1=0.025, y1=0.03, x2=0.025, y2=0.0, gapx=1)
# print(mi)

def calculate_mutual_inductance(x1, y1, x2, y2, gapx):
    """
    used to get machine parameters from machine.xlsx file
    path ='C:\\Users\\xiang\\Documents\\XINAO\\AI\\SEPY\\data\\EXL50U\\machine.xlsx'

    input:
        option name         data_type    size    unit    meaning
        1      x1           float        N               x-coordinate of a point in the coil L1
        2      y1           float        N               y-coordinate of a point in the coil L1
        3      x2           float        1               x-coordinate of a point in the coil L2
        4      y2           float        1               y-coordinate of a point in the coil L2
        5      gapx


    output:
        option  name                data_type   size    unit    meaning
        1       mutual_inductance   folat       N       H       mutual inductance between two current elements

    reference:
        [1]皇甫国庆.两圆线圈间互感及耦合系数讨论[J].渭南师范学院学报,2015:26-31.

    Developed by xiangguchn 2023/11/28
    If you have some question, don't hesitate to ask xiangguchn@163.com for help.
    """
    cmu = np.pi * 2.0e-7

    R1 = np.sqrt((x1 + x2) ** 2 + (y1 - y2) ** 2)
    m = 4 * x2 * x1 / R1**2
    
    myk = special.ellipk(m)
    mye = special.ellipe(m)
    mutual_inductance = cmu * R1 * (2 * (myk - mye) - m * myk)

    return mutual_inductance


if __name__ == "__main__":
    mi = calculate_mutual_inductance(x1=0.025, y1=0.03, x2=0.025, y2=0.0, gapx=1)
    print(mi)
    test_calculate_mutual_inductance(x1=0.025, y1=0.03, x2=0.025, y2=0.0, gapx=1)
